%% fig2
clc; clear; close all;

%% —— 参数设置 —— 
g11        = 1.4;
g12        = -2;
g21        = 1.4;
g22        = 1.6;

dt         = 0.01;      % 时间步长
t_total    = 300;       % 总积分时间
renormStep = 50;        % 每隔 renormStep 步做一次正交化

% 扫描范围：k ∈ [−8,8], x₃₀ ∈ [−3,3]
nk       = 40;
n30      = 40;
k_list   = linspace(-8,  8, nk);
x30_list = linspace(-3,  3, n30);
[K, X30] = meshgrid(k_list, x30_list);

% 预分配前两条 LE 曲面
LE1 = zeros(size(K));
LE2 = zeros(size(K));

%% —— 并行计算 LE₁ 和 LE₂ —— 
parfor idx = 1:numel(K)
    k   = K(idx);
    x30 = X30(idx);
    x0  = [0.01; 0; x30];  % 初始状态 (x1₀,x2₀,x3₀)

    % 系统方程与雅可比句柄
    f = @(x)      MBNHNN(x, k, g11, g12, g21, g22);
    J = @(x)      Jaco(  x, k, g11, g12, g21, g22);

    % Benettin 方法计算前三个 LEs
    Ls       = Benettin(f, J, x0, dt, t_total, renormStep);
    LE1(idx) = Ls(1);
    LE2(idx) = Ls(2);
end

%% —— 绘图 —— 
figure; hold on;

% LE1 曲面
h1 = surf(K, X30, LE1, 'EdgeColor','none', 'FaceAlpha',0.8);
% LE2 曲面
h2 = surf(K, X30, LE2, 'EdgeColor','none', 'FaceAlpha',0.8);
% LE1 投影
h3 = surf(K, X30, -ones(size(K)), LE1, ...
    'EdgeColor','none', 'FaceAlpha',1);

% 自定义 colormap：白→蓝→黑→红→黄→白
pos  = [0, 0.25, 0.5, 0.75, 0.9, 1];
cols = [1 1 1; 0 0 1; 0 0 0; 1 0 0; 1 1 0; 1 1 1];
xi   = linspace(0,1,256)';
cmap = interp1(pos, cols, xi, 'linear');
colormap(cmap);

% 颜色范围 ±LE1 最大绝对值
M = max(abs(LE1(:)));
caxis([-M, M]);
colorbar;

view(-45,30);
xlabel('k',       'FontSize',12);
ylabel('x_{3,0}', 'FontSize',12);
zlabel('LEs','FontSize',12);

legend([h1 h2 h3], {'LE_1','LE_2','LE_1 投影'}, 'Location','best');
grid on;
axis tight;

